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Abstract. We have performed Smoothed Particle Hydrodynamics (SPH) simulations to study the time evolution of one and two 
protoplanets embedded in a protoplanetary accretion disc. We investigate accretion and migration rates of a single protoplanet 
depending on several parameters of the protoplanetary disc, mainly viscosity and scale height. Additionally, we consider the 
influence of a second protoplanet in a long time simulation and examine the migration of the two planets in the disc, especially 
the growth of eccentricity and chaotic behaviour. One aim of this work is to establish the feasibility of SPH for such calcula- 
tions considering that usually only grid-based methods are adopted. To resolve shocks and to prevent particle penetration, we 
introduce a new approach for an artificial viscosity, which consists of an additional artificial bulk viscosity term in the SPH- 
representation of the Navier-Stokes equation. This allows accurate treatment of the physical kinematic viscosity to describe the 
shear, without the use of artificial shear viscosity. 
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1. Introduction 



Since the discovery of the fi rst extrasolar planet in t he mid- 
nineties of the last century by M ayor & Ouelo j dl995h . the in- 
terest in planet formation research has grown rapidly. By now 
about 120 extrasolar planets are known. For comprehensive 
reviews on extra sola r planets and their detection we refer to 
IPerrvmanl d2000h and lMarcv et alJll2000h . These newly discov- 
ered planetary systems partially feature attributes that differ 
enormously from those of our solar system. In particular, or- 
bital semi-major axes of down to 2.25 x 10~ 3 AU, masses in 
the range of 0.12 - 13.75 Mj up , and large eccentricities of up 
to 0.927 have been found, in contrast to the planets in the solar 
system which orbit the sun at larger distances on almost circu- 
lar orbits with rather lower masses. Therefore, nowadays planet 
formation theory not only has to explain the formation of our 
solar system, but the formation of planetary systems with these 
substantially different features as well. 

It is generally believed that planet formation is part of 
star formation. Stars are formed by gravitational collapse of 
a molecular cloud. During this process, angular momentum 
conservation leads to the formation of a circumstellar disc 
around a central object, the protostar. The material in the 
disc slowly accretes onto the protostar in the centre until the 
protostar becomes a star. Circumstellar discs composed of 
dust and gas have been fo und around young T Tauri stars 
jBeckwifh & Sargeiill [l996). Star formation theory suggests a 



typical lifetime of such circumstellar discs in the range of 10 6 
to 10 7 y r, which determines the time frame for planet forma- 
tion (see Kallenbac h et alJ J200fj|) and references therein). It is 
assumed that planets form in these circumstellar (or protoplan- 
etary) discs. Two mechanisms are proposed, either fast forma- 
tion by direct fragmentation due to gravitational instabilities 
( Mave r et al]200alBosj200lll2003l) . or slow formation during 
a process of several steps: at first the dust in the disc grows via a 
so-called hit-and-stick mechanism and settles in the mid-plan e 
of the disc, where a dense layer forms ( Beckwith et al. 2000). 
There the solid material is able to combine into larger bod- 
ies, which eventually form planetesimals. Finally, planetary- 
sized objects d evelop out of these planetesimals by collisions 
JOhtsuki et alJ 120021) . The latter, slowe r scenario lea d s to a 
longer formation time. Calculations by Polla ck etal] {1 996) 
suggest a formation time of 1 to 10 million years for Jupiter 
and Saturn, and of 2 to 16 million years for Uranus, respec- 
tively. 



This scenario leads to essentially circular orbits of the plan- 
ets as they form in a keplerian disc. Due to the lack of material 
in the inner regions of the disc, more massive planets can only 
be created at larger distances to the central object. However, 
while these properties fit well to our solar system where the ec- 
centricities are rather small and the massive planets are located 
in the outer regions, the model cannot explain the appearance 
of highly massive planets on significantly eccentric orbits at 
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orbital radii down to 0.04 AU, as found among the extrasolar 
planets. 

A possible solution to this problem is planet migration. It is 
generally assumed that these planets have formed further out 
in the disc and have migrated inwards later. This migration 
has its cause in tidal in teractions of the planet with the disc 
jGoldreich & Tremainelll980l) . The planet exerts tidal torques 
on the disc which lead to the formation of spiral density wave 
perturbations. Simultaneously, since the disc loses its axisym- 
metry, a net torque acts o n the planet itse lf, which starts mi- 
grating in radial direction (Lin et alfc oOO). The migration rate 
depends on different factors, such as the mass of the planet, 
the density distribution in the disc and the viscosity of the disc 
material. 

Another problem lies in the high masses of the discovered 
planets. The tidal interaction of the planet can eventually lead 
to the opening of a gap at the planet's orbital region in the 
disc, which has a severe im pact on the accretion rate. Un til the 
numerical calculations bv lArtvmowicz & Lubowl II1996). who 
modelled the accretion on circumbinaries, and t he simulation s 
of a iovian p rotoplanet embedded in a disc by 
iBrvden et alJ Jl999h . it has been assumed that accretion stops 
after gap formation. These simulations however yield a final 
non-vanishing accretion rate which depends on the viscosity 
in the disc. Since then, furt her work has been done on this by 
various other authors JLubow et alJI 19991 IBrvden et alJ feoOO: 
Lin et al.l2000llTerauem et alJ2000HNelson et al.l200ol) . 

While the numerical method known as Smoothed Particle 
Hydrodynamics (SPH) has been used for one of the first sim- 
ulations t hat examined the sy stem of an embedded protoplanet 
in a disc ( Seki va et al.ll9 87). later simulations were nearly ex- 
clusively performed with grid-based methods, mainly because 
computer resources have not been available for other high res- 
olution simulations. This has changed only recently. In this pa- 
per we present a feasibility study for the simulation of planet- 
disc interactions with SPH, applying a new approach for the 
artificial viscosity, and we show the results of simulations with 
one and two embedded protoplanets in a viscous accretion disc. 
Very recently, also iLufkin etal have used SPH for in- 

vestigating planet-disc interaction. However, they focus on the 
scenario where planets form due to selfgravitating instabilities, 
while we consider non-selfgravitating discs in two dimensions 
(r-(p). 

The outline of the paper is as follows. In the next sec- 
tion we present the physical model of the system, including 
an overview on gap formation and planet migration theory. 
Numerical issues will be discussed in Sect. [3] where espe- 
cially the approach for a new artificial bulk viscosity will be 
described. The results of several simulations will be given in 
|4] followed by a discussion in Sect. [5] Finally, we draw our 
conclusions in Sect. [5] 

2. Physical model 

We only consider a non-selfgravitating, thin accretion disc 
model for the protoplanetary disc. The disc is located in the 
z — plane, rotating around the z-axis. The vertical thickness 
H of the disc is assumed to be small in comparison to the radial 



extension r, H/r «: 1. Therefore it is possible to vertically inte- 
grate the hydrodynamic equations and use vertically averaged 
quantities only, neglecting the vertical motion in the disc. The 
protostar in the centre of the disc always is assumed to have a 
mass of M c = 1M©. 

2.1. Basic equations 

The equations describing a viscid flow of gas in the disc are 
the continuity equation and the Navier-Stokes equation. In the 
Lagrangian scheme the former reads 



^ + EVv = 0, 
dt 



(1) 



where S is the surface density and v the velocity. The surface 
density is defined by 



-r 



Qdz, 



(2) 



where g denotes the density. The Navier-Stokes equation in 
Lagrangian formulation reads 



^ = _I V p +/+ IvT, 

dt S £ 



(3) 



where p is the vertically integrated pressure, and / denotes a 
specific external force, i.e., the gravitational acceleration due 
to the star and the N p planets with masses m, and locations r,- 



/ = ~GM Q 



- r.V Zj 



(r - r c ) 



- V Gnti- 



(4) 



The last term in the equation of motion describes the viscous 
forces and contains the viscous stress tensor which can be writ- 
ten as 

<9v. 



.dv p dv a 2 <9i 



dx a dxn 3 dx 



8Vy 
8Xy 



(5) 
(6) 



(Greek indices denote spatial coordinates, and Einstein's sum- 
ming convention holds throughout), where the coefficients of 
shear viscosity, 77, and bulk viscosity, are positive and do not 
depend on the velocity. 

2.2. Equation of state 

For computational simplicity we follow the usual approach and 
use an isothermal equation of state, where the vertically inte- 
grated pressure p is related to the surface density S by 



P = <s. 

Here the local isothermal sound speed is given by 



H 

C s = — Vkep, 



(7) 



(8) 



where v^p denotes the Keplerian velocity of the unperturbed 
disc. The scale height of the disc, which is defined as the ra- 
tio of the vertical height to the radial distance, H/r, equals the 
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inverse Mach number Ai of the flow in the disc and is taken 
as an input parameter for our disc model. By choosing such an 
equation of state we suppose that the disc stays geometrically 
thin and emits all thermal energy which results from tidal and 
viscous dissipation locally. 

2.3. Viscosity 

Viscous processes are of great importance in the theory of 
accretion discs. In order to move radially inwards, the disc 
material has to lose its angular momentum. Therefore a per- 
manent flow of angular momentum radially outwards and of 
mass radially inwards takes place in the disc. Often, the- 
ories of protostell ar discs assu me viscous accretions discs, 
see e.g. Papaloizou et al We follow the prescription 

by IShakura (felsunvaevl i 19731) and assume an effective a- 
viscosity, where the kinematic viscosity v in the disc is given 
by 



v = ac s H. 



(9) 



The shear viscosity rj is related to the kinematic viscosity by 
r] = VL. Typical a-values for protostellar discs are in the range 
of 1(T 2 to 1(T 3 . 



2.4. Gap formation and planet migration 

Using perturbation theory, the tidal interaction of the planet 
and the accretion disc can be analyzed analytically if the plan- 
etary ^jnas^jsjnu£h_smaUer thanthe_mass of the central object, 
see Goldr eich & Tremaind dl979i Il980h for details. The planet 
exerts torques on the disc material (and vice versa), which 
can lead to gap opening at the planet's orbital region if these 
tidal torques exceed the internal viscous torques in the disc 
(Papaloizou & Lin 1984). The faster (slower) moving disc ma- 
terial inside (outside) the planet's orbit exerts a positive (neg- 
ative) torque on the planet. The planet will begin to migrate 
radially as soon as these torques do not cancel out. There are 
two different timescales for migration depending on whether 
a gap has formed or n ot, which a re called type I and type II 
migration respectively (Ward 1997). 

Type I migration. No gap formation occurs as long as the 
response of the disc to the tidal perturbation stays linear. 
The radial migratio n velocity of the planet is then given by 
llNelson et al l2000l) 

d7~-^(£)(^)(^) 3 ' (10) 

where q denotes the mass ratio of the planet to the central star, 
M c the mass of the central star, and Q.^ ep the Keplerian angu- 
lar velocity. The factor c\ accounts for the imbalance of the 
torques between the two regions of the disc inside and outside 
the planet's orbit. Hence, the radial migration velocity will be 
larger for a more massive planet as long as the response of the 
disc stays linear. 

As soon as the planet has grown to a sufficient high mass, 
the response of the disc to the perturbation becomes non-linear, 
and gap formation occurs. 



Type II migration. In order to form a gap the tidal torques 
exerted by the planet have to exceed the internal viscous 
torques in the disc. The conditio ns for gap formati o n hav e 
been widely stud ie d in d etail b y |Pa^aMzou & Lirj (H'984 ), 
| Lin & Papaloizou! dl985l 1 19931) . iBrvdenetalJ i 19991) and 
Nel sonetalT(i200ol) . Here we only present the different trunca- 
tion conditions and scenarios for type II migration. A sufficient 
condition for viscous truncation is 



40v 
q > — j . 



(11) 



where to is the angular velocity of the planet and r 
its orbital radius. For discs w ith a very low viscosity 
iKorvcanskv & Papaloizou! \ 19961) suggest the thermal trunca- 
tion condition 



q> 



(12) 



During type II migration, the migration of the planet is driven 
by the viscous evolution of the disc. The migration rate in this 
case is given by the radi al drift velocity of the g as in the disc 
due to viscous processes (Papaloizou & Lin 1984) 



dr 
It 



3v 
2r" 



(13) 



This relation holds as long as the mass of the planet is lower 
than the mass of the disc material with which it interacts. 
Otherwise the inertia of the planet gets important and slows 
do wn the orbital evolut ion. This has been studied in d etail 
bvlSver&Clarkel(ll995l) and llvanov etaDjl999fc . lNelson et alJ 
( 2000h find a migration rate based on a model bv llvanov et alJ 
( 1999) which reads 



df 2 v I p 

where m p is the mass of the planet. 



(14) 



3. Numerical issues 

We apply Smoothed Particle Hydrodynamics (SPH) for 
our simulations. SPH is a grid-free Lagrangian particle 
method for solving the system of hydrodynamic equa- 
tions for compressible and v isco us fluids first introduced by 
iGingold & Monaghan] i ll 977b and iLucvl d 1977b . SPH is espe- 
cially suited for boundary-free problems with high density con- 
trasts. Rather than being solved on a grid, the equations are 
solved at the positions of the so-called particles, each of which 
represents a volume of the fluid with its physical quantities like 
mass, density, temperature, etc. and moves with the flow ac- 
cording to the equations of motio n. For a more det ailed de- 
scription see, e.g.. lBenz] (fl990) and Monaghanl i ll 9921) . 



3.1. The code 

Our SPH-code is fully parallelized to make optimal use of to- 
day's supercomputers which is necessary to achieve the re- 
quired resolution and accuracy. The parallel implementation is 
done with the so called ParaSPH Library, a generalized set of 
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routines to simplify and optimize the development of parallel 
particle codes (Bu beck et alJIl 9991) . This library allows a clear 
separation between the implementation of the physical problem 
to be simulated and the parallelization. A programmer using the 
library works with parallel arrays containing the physical quan- 
tities of the SPH-particles and iterators to step through the par- 
ticles and their neighbours. Domain decomposition, load bal- 
ancing, nearest neighbour search and communication is hidden 
in the library. ParaSPH uses MPI (Message Passing Interface) 
for the communication and works on every parallel architec- 
ture providing an MPI library. On a Cray T3E one can achieve 
a parallel speedup of about 120 using 256 nodes for a prob- 
lem with 360000 particles. On the Beowulf-Cluster where we 
have performed the majority of our simulations - a Pentium III 
based Linux-Cluster with a Myrinet communication network 
located at the University of Tubingen - the parallel speedup is 
about 60 on 128 processors. 

3.2. Physical and artificial viscosity, XSPH 

To model the physi cal shear viscosity w e use the implementa- 
tion introduced by Fleb be et alJ {l994) and solve the Navier- 
Stokes equation including the entire viscous stress tensor (|6j, 
in contrast to the usual approach of a scalar artificial viscosity 
term. For the physical shear flow, we assume a constant kine- 
matic viscosity v and a vanishing bulk viscosity £ = 0. Then, 
the acceleration due to the shear stress tensor cr a n reads 



dvv 
df 



shear 



1 5(vEcr ff/ j) 

2 dxn 



(15) 



Applying the smoothing discretization scheme to the derivative 
of the stress tensor yields the SPH representation 



dVa, 

df 



shear 



j 



dx B 



with the particle form of the shear tensor 

C" aBi = V a Ri + Vfl, 



2 

3* 



afi ** 7 yyi j 



where V a pi is the particle representation of the velocity gradient 

dva/dxp 



VaBi - 



dxn 



-(Vaj ~ Vai) 



dxa ' 



The indices ; and j denote quantities that are assigned to SPH 
particles with positions r, and r,, respectively; ntj is the particle 
mass, and W is the usual SPH kernel function. Replacing sur- 
face density S by (volume) density g, this SPH-representation 
of the shear viscosity also holds in 3D. 

In order to prevent the particles from mutual penetration, 
we add an additional artificial bulk viscosity to the viscous 
stress tensor. The bulk part of the viscous stress tensor is given 
by £Vv and leads to an additional acceleration of particle i ac- 
cording to 



dv, 

df 



bulk 



= 2 m J& 



(VV),- + (VV); 



The artificial bulk viscosity coefficient for the interaction of 
particles i and j is given by 



/,.. = ] -A ; <Vv ),;!,,, forvyry<0 
1 else 



(20) 



where hi (hj) is the smoothing length of particle i (j), Q,j is an 
abbreviation for the average (Q-, + Qj)/2 for any quantity Q, 
and the notation Vyfy = (v, - vj) ■ (r, - rj) has been used. Only 
approaching particles are contributing to the artificial viscosity. 
The factor / is of order unity and set to 0.5 in our simulations. 
The divergence of the velocity is calculated according to 



(Vv) ' = ^Z m > (v '- v > )VW '> 



(21) 



To avoid mutual particle penetration, we additionally use 
the XSPH-algorithm. The particle s are moved at av eraged ve- 
locity of their interaction partners cMo naghanl 19891) 



An 
df 



= V; + X 



z 



2m i 



(Vi-VjWij, 



(22) 



where x is in the range between and 1 . The use of this aver- 
aged velocity keeps the particles in their relative local order. 

In appendix|A]we demonstrate the low intrinsic numerical 
shear stress of this new artificial bulk viscosity approach. 

3.3. Migration 

To study the migration of the protoplanet in detail, the equation 
of motion of the protoplanet is integrated. The planet's acceler- 
ation due to the gravitational forces reads 



(16) ^1 



df I- _ r I 



r, 



N r - 
^ Gm, — 

<' =1 Ik - r,f + 



\3/2 



(23) 



(24) 



where M c is the mass of the central object, r c its position, Af 

(17) the total number of SPH particles in the simulation with posi- 
tions r,-, and G the gravitational constant. Only particles which 
are not inside the Roche lobe of the planet contribute to the 
sum. Additionally, the gravitational forces of the particles on 
the planet are softened by the use of 6 2 , which is s mall. The 

(18) Roche lobe is approximated using the formula by lEggletonl 
ill 9831) 

0.49^ 2/3 
rroche " 0.6 9 2/3 + ln(l+? 1/3 )'' P ' 

In the case of more than one planet, the additional gravitational 
forces of the other objects are added to Eq. J23I . nevertheless 
this approximation formula for the Roche lobe is used. 

3.4. Initial conditions 

The initial setup of our standard model places a jovian planet 
with mass m p = 1 Mj up at r p = 5.2 AU and a central object with 
(jo,) mass M c - 1 M©, which gives a mass ratio of q — 9.55 x 10~ 4 . 
The surface density 2 falls off as r~ 3/4 , and the particles with 
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equal masses are distributed according to this density between 
'"min = 1 AU and r max = 10 AU. The mass of the disc is set 
to 1CT 2 Mq and the initial particle velocities are v r = and 
v v - v kep - >/GM c /ri as in the unperturbed case. The typical 
number of particles at the beginning of the simulation is 3 x 10 5 . 
The standard model uses a Mach number of M — 20, that is a 
scale height H/r = 0.05, and the kinematic viscosity parameter 
vo is set to 10 15 cm 2 /s, which corresponds to an a value of 
about 4 x 10~ 3 at the orbital radius of the planet. 

As the density falls off with radius, fewer particles are lo- 
cated in the outer regions of the disc than near the central ob- 
ject. Hence, a constant number of interaction partners is more 
suited to the problem than a fixed smoothing length. Unless 
otherwise mentioned the number of interaction partners per 
particle is set to 100. As mentioned above, we use the XSPH- 
method to avoid mutual particle penetration and we set x = 0.5 
in all our simulations. 



4. Results 

In this section we present the results of the simulations with 
various input parameters. First we consider the case of a single 
planet in the disc. Afterwards, the case of two planets and their 
mutual interactions is investigated. 

4.1. One planet 
4.1 .1 . Gap formation 

In order to study the formation of a gap, we have performed 
several simulations of a disc with an embedded protoplanet. 
The disc scale height and the kinematic viscosity have been 
used as input parameters and have been varied accordingly. 
Here, the gravitational forces of the particles on the planet have 
not been considered, keeping it on a circular orbit around the 
central star. Moreover, the planet was assumed to be able to 
accrete material without growing in mass, its Roche lobe has 
been kept constant throughout the simulations. 

Just after the start of the simulation the planet disturbs the 
density distribution in the disc effectively, which leads to the 
typical spiral density wave pattern in the disc which can be 
seen in Fig. ^ where a greyscale plot of the surface density 
is presented for the system after 10 orbits of the protoplanet. 
In the case of this reference model, gap formation begins im- 
mediately after the start of the simulation run. After ten orbits 
of the planet the decrease of the density at the protoplanet's 
orbital region is already clearly visible. The evolution of the 
azimuthally averaged surface density during the gap formation 
process is shown in Fig. [2] After one hundred orbits of the pro- 
toplanet, the decrease of density at the orbital region is more 
than two orders in magnitude. 

The radial size of the gap depends only slightly on the vis- 
cosity, as can be seen in Fig. [3] where the surface density is 
plotted for three simulation runs with different kinematic vis- 
cosity: v = vo = 10 15 cm 2 /s (reference model), v = 2vo, and 
v = 4vo- A higher viscosity leads to a faster radial motion of 
the gas in the disc, and thus to more loss at the inner bound- 
ary, which results in a slightly lower density profile after the 



3. 

us 




nJ2 



[rad] 



2* 



Fig. 1. Greyscale plot of the surface density in the disc in the 
case of the reference model after 10 orbits of the protoplanet 
in the disc. The gap clearing has already started, and the spiral 
density wave pattern is visible throughout the disc. The cir- 
cle represents the Roche lobe of the protoplanet according to 
Eggleton's formula. 



600 



CD 
CD 



after 1 orbit 
after 1 orbits 
after 20 orbits 
after 50 orbits 
after 1 00 orbits 



500 



400 



300 



200 



100 



4 6 8 10 

Distance to central object [AU] 

Fig. 2. Evolution of the azimuthally averaged surface density 
for five different times in the case of the reference model. The 
planet has a fixed orbit at 5.2 AU. 




same n umber of orbits. A similar behaviour was found bv lKlevI 
Because in Fig. [3] the simulation time is rather short 
compared with the viscous timescales, the gap structure looks 
very similar in each case. The differences may become more 
pronounced for longer simulation runs. Unlike the viscosity, 
the scale height has a severe impact on the depth of the gap as 
can be seen in Fig.0] There, the surface density of the disc af- 
ter 100 orbits of the protoplanet is plotted for the scale heights 
H/r = 0.1, H/r = 0.05, and H/r = 0.025. The smaller the 
scale height, the lower the temperature and the deeper the gap, 
as can be seen clearly. Additionally, the spiral density wave pat- 
tern gets tighter as the sound speed is lower for smaller scale 
heights. 
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2 4 6 8 10 12 14 

Distance to central object [AU] 

Fig. 3. Influence of the kinematic viscosity parameter on the 
gap size. The azimuthally averaged surface density is plotted 
for three different viscosities in units of vrj = 10 15 cm 2 /s after 
100 orbits of the protoplanet. The planet has a fixed orbit with 
a radial distance of 5.2 AU. 



4.1.2. Migration and accretion 

In order to study the migration of the protoplanet we consider 
again the standard reference model of a jovian protoplanet in 
a circular orbit with radius 5.2 AU. To achieve a steady-state 
particle distribution we proceed as in the last section and let the 
system evolve until the gap is nearly fully developed before we 
allow the gravitational force of the particles to act on the proto- 
planet. Then the simulation is restarted and the trajectory of the 
planet is integrated as well. For the reference model we expect 
pure type II migration on a viscous timescale. For this simula- 
tion we used nearly 360 000 particles to model a disc with mass 
3.33 x 10~ 3 M and r m = 1 AU, r out = 13 AU. Particles that 
get closer to the protoplanet than 50 % of its Roche radius are 
taken out of the simulation and assumed to be accreted by the 
protoplanet. In one model, the mass of the accreted particles is 
added to the mass of the planet, in a second, the mass of the 
planet is fixed for all times. The total simulation time is about 
30000 years. 

The evolution of the radial distances of the planets is shown 
in Fig. |5] Initially the migration rate is very high and about the 
same for both the planet with fixed mass and the planet with 
increasing mass. After 3 000 years the radial distances to the 
central object differ distinctly, with the more massive planet 
migrating slightly faste r. Similar behaviou r was found in nu- 
merical simulations by Nel son et alJ (120001) . The planet with 
a fixed mass ends up at an orbital distance of 4.5 AU after 
20 000 years, which would correspond to a steady-state migra- 
tion rate of 3.5 x 10~ 5 AU/yr. The planet with increased mass 
has a radial distance of 4.27 AU after 30 900 years, or a steady- 
state migration rate of 1.4 x 10~ 4 AU/yr. As can be seen in 
Fig. |3 the migration rates are not constant but decrease con- 
siderably. This is caused primarily by the mass loss of the disc, 
leading to less gravitational interaction between the planet and 
the disc. 




*12 Tt. 2tl 



H/r=0.Q25 




- " 



Fig. 4. Influence of the scale height of the disc. The greyscale 
plot shows the surface density in the disc after 100 orbits of the 
protoplanet for three different scale heights. 



The accretion rates of the planets for both models are pre- 
sented in Fig. |6] The initial accretion rate of about 2.5 x 
10" 4 Mj U p/yr drops rapidly in the first 3 000 years to a value of 
5 X 10~ 5 Mj U p/yr for both the planet with constant mass and the 
planet with increasing mass. At that time, the inner part of the 
accretion disc has vanished, as the material was either accreted 
by the protoplanet or by the central object. After the loss of the 
inner disc, the accretion rate decreases more slowly. Between 
3 000 and 7 500 years of simulation time, the planet with in- 
creasing mass has a slightly lower accretion rate. The situation 
changes after 10 000 years, when the planet with the fixed mass 
accretes less. At the end of the simulation, both planets have 
about the same accretion rate of 10~ 5 Mj llp /yr, although the 
planet with mass accumulation has already doubled its mass. 
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Fig. 6. Evolution of the accretion rate on the planet for a planet 
with fixed mass and a planet with increasing mass. The planet 
is originally located at a radial distance of 5.2 AU, and migrates 
inwards as shown in Fig. EI 



However, during that time the constant loss of material through 
the open boundaries and the accretion on the planet lead to a 
decrease in disc mass of nearly one magnitude, which has to 
be considered for comparisons to grid-based simulations with 
closed boundaries. 



4.2. Two planets 

As it is generally assumed that the mutual interaction of grav- 
itating objects in an accretion disc can finally lead to eccentric 
orbits, we have included another protoplanet into our simula- 
tions. 



4.2.1. Gap formation 

For the simulation with two protoplanets the setup was chang- 
ed as follows: one protoplanet is located at a distance of 5.2 AU 
from the central gravitating object, the other one at twice that. 
The kinematic viscosity parameter is again set to 10 15 cm 2 /s, 
and the disc scale height is 0.05. Now the outer radius of the 
disc is 20 AU, and the inner boundary is set to 2AU. The num- 
ber of particles in this run was 320 000 and the disc mass was 
set to 10~ 2 M0. The evolution of the density in the disc is 
shown in Fig. As in the case of a single jovian planet, gap 
formation starts immediately after the beginning of the simula- 
tion. This happens for each planet individually, until the parti- 
cles between the inner and outer planet have been accreted by 
the inner planet, and a common gap has formed. The innermost 
region of the disc has already vanished after 200 orbits of the 
inner planet. 

4.2.2. Migration and accretion 

After 200 orbits of the inner planet, the simulation is restarted 
and the trajectories of the two planets are integrated to exam- 
ine their migration. The evolution of the semi-major axis of the 
two planets is shown in Fig.|S]for the next 30 000 years. During 
the first 5 000 years the inner planet mainly stays on its orbital 
radius around the central star because it is not surrounded any- 
more by disc material (see Fig. [7} which could exert torques. 
On the other hand the outer planet migrates inwards due to the 
action of the outer disc, resulting in a radial approach of the two 
planets. As the orbital distance of the two planets decreases, 
they capture each other in a 2:1 mean motion resonance at ap- 
proximately t — 7 500 years. Upon capture their eccentricities 
slowly increase, as can be seen in Fig. [9] By the end of the 
simulation, the eccentricity of the outer planet has grown to 0.1 
and the eccentricity of the inner planet to 0.18. In the end, the 
planets have reached masses of 1 .34 (inner planet) and 2.9 Mj up 
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Fig. 9. Evolution of the eccentricities of the two protoplanets 
during the migration, same simulation as presented in Fig. [8] 



(outer planet), respectively. The mass of the disc changed from 
1(T 2 M at the start of the simulation to 5.64 x 10 3 M at the 
end. After all, only 20 % of the mass lost by the disc have been 
accreted by the protoplanets. The resonant capture and overall 
evolution of this system is nearly identical to that described in 
lKlevl(l2000l) . 

5. Discussion 

In the case of a single jovian planet we obtain the same influ- 
ence of the disc scale height on the gap formation and on the 
structure of the spiral wa ve pattern as found already by se veral 
other authors (see, e.g., iBrvden et all 1 19991 lKlevlll999l) . and 
as postulated by the theory. For the parameters of the reference 
model, Eq. (1131 yields a migration rate of 4.13xl0~ 5 AU/yrfor 
pure type II migration, for a planet located at a radial distance 
of 5.2 AU. We find migration rates of the same order of magni- 
tude by averaging over the complete simulation time. The mi- 



gration rate for this simulation agrees well with migration rates 
found bv lD' Angelo et all d2003l) . who examined the migration 
and growth of one protoplanet in a disc for various initial planet 
masses, using a grid-based method and closed boundary condi- 
tions. 

Interestingly, the accretion rate for a planet with a fixed 
mass and a planet with increasing mass do not differ con- 
siderably at the end of a simulation time of 30000 years. 
Presumably due to the loss of disc material at the open bound- 
aries, the resolution is too coarse to resolve the difference 
correctly. After a simulation time of 15 000 years, the accre- 
tion rate onto the planets is about 2.5 x 10~ 5 Mj up /yr for 
our reference model with a kinematic viscosity coefficient of 
10 15 cm 2 /s. This val ue agrees well with t he results from simu- 
lations performed by Lubo w et al 1 dl999l) . who used the ZEUS 
hydrodynamics code to simulate disc accretion onto high-mass 
planets. They find an accretion rate of 2.13 X 10~ 5 Mi un /yr , 
which is half the valu e that was published by iKlevI (1999). 
However, iLubow et alJ dl999h use different boundary condi- 
tions, namely reflecting closed boundaries at the inner edge of 
the disc, and open boundaries at the ou ter edge. The results of 
the simulations by Brvd en et alJ (11999) provide accretion rates 
in the same range, although they use an initial density profile 
which is constant throughout the disc. Extensive simulations 
with three different grid codes, nirvana, rh2d, and f argo, by 
iNelson et alJ |2000) provide comparable accretion rates for a 
jovian planet. 

For a single planet in the accretion disc the orbit remains 
circular, although the planet migrates towards the central star. 
This suggests that the gravitational interaction between the 
planet and the protoplanetary disc cannot excite a change in 
the ecce ntricity of the orbit of t he planet, at least for low mass 
planets ( Papaloizo u et al 1 12001 ). Therefore, in systems with a 
single extra-solar planet mostly circular orbits should be ex- 
pected. This clearly appears to be in contradiction to the ob- 
servations, and is why different mechanisms, e.g. resonances, 
for the evolution of high eccentricities in single-planet systems 
have to be assumed. 

We have found that a number of 360000 SPH-particles, 
which has been used in the simulation of migration and planet 
accretion, is sufficient in 2D to resolve an expected accretion 
rate of several 10 Mj up /yr. However, the use of an artificial 
bulk viscosity is essential, as extensive test calculations have 
shown. Without the use of this additional viscosity, the particle 
orbits penetrate each other, density shock waves cannot be re- 
solved, and even the formation of the gap may be suppressed. 
Note that in contra st to the normally used a rtificial a- viscosity 
approach (see Monaghan & Gingold ( 1983) ) which le ads to an 
effective physical shear viscosity (Nel son et al.lll998l) . this ar- 
tificial bulk viscosity has very little impact on the shear kine- 
matic viscosity in the disc, as can be seen in appendixlAl 

The number of SPH-particles was not high enough though 
to study the flow arou nd the protoplanet wit hin its Hill ra- 
dius, as was done by D'Ange loetall ^002) using nested- 
grid techniques. This would have made it possible to exam- 
ine the accretion proce ss onto the protoplanet more precisely. 
ID' Angelo et al.l d2002 ) find the formation of a so-called proto- 
jovian disc around the protoplanet which dominates the accre- 
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tion process onto the planet. However, in our simulations the 
spatial resolution was too low to resolve the formation of such 
a disc. 

Generally, SPH seems to be slightly more computationally 
expensive for this kind of application than comparable simula- 
tions with grid-based schemes. This is mostly due to the intrin- 
sic noise of the particle method that requires a comparatively 
high particle number to achieve a similar spatial resolution. The 
advantage of this numerical method is its Lagrangian nature, 
which allows an easy treatment of open boundaries and com- 
plex geometries, like multiplanetary systems. 

In order to study the evolution of a protoplanetary system 
with several planets, we included a second protoplanet in the 
simulation. After the first 200 orbits of the inner planet, the 
inner disc has already nearly vanished leading to a lower ac- 
cretion and migration rate of the inner planet in comparison to 
the outer planet. Due to different migration speeds of the two 
planets they approach each other radially and are eventually 
captured into a 2:1 mean motion resonance. Upon capture, the 
dynamical interaction between the two planets increases and 
their eccentricities begin to grow considerably. This resonant 
behaviour has been fou nd and investigated numerically already 
by several authors (e.g. [ Kiev 200o l]Nelson & Papaloizour2 002; 
Lee & Peale 2002; Kiev et al. 2004). Depending on the subse- 
quent evolution (distance of migration and damping), the sys- 
tem might become unstable eventually. 

Thus, the gravitational interaction between protoplanets in 
the disc can explain the appearance of resonant extra-solar 
planets with high eccentricities as found for example in GJ 876, 
HD 82943, or 55 Cnc. 

6. Conclusion 

With the presented calculations we have modelled the dynam- 
ical evolution of a protoplanet embedded in a disc using the 
Lagrangian particle method SPH. The SPH-representation of 
the shear viscosity allows accurate treatment of the viscosity of 
the system, and the use of a newly introduced artificial bulk vis- 
cosity term in the Navier-Stokes equation inhibits mutual par- 
ticle penetration. Our investigations comprised the evolution of 
a single Jovian-like protoplanet in a protoplanetary accretion 
disc, including the accretion of gas onto the planet and the mi- 
gration of the planet through the disc. Moreover, we examined 
the effects on migration if two planets are located in the disc. 

According to our simulations, gap formation does not in- 
hibit accretion onto the protoplanet. Our simulations indicate 
that a protoplanet can grow up to several Jupiter masses be- 
fore it has migrated to the central star. The simulations with 
two protoplanets show a possible mechanism for the formation 
of a two-planet system with a lower-mass planet near the cen- 
tral star and a more massive planet with a larger semi-major 
axis. The results of our SPH simulations compare favourably 
with those of grid codes. Although the computational costs are 
rather high, and therefore SPH seems to be less effective to 
model this problem, it is well suited to verify simulation results 
achieved with other, grid-based numerical schemes because of 
its completely different numerical approach. Moreover, it may 
be useful for 3D-simulations or self-gravitating flows. 
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Appendix A: Numerical shear 

To measure the intrinsic numerical shear viscosity of the new 
bulk viscosity approach, we have performed several test simu- 
lations of the viscously spreading ring, a model representing an 
idealised accretion disc orbiting a point mass. Because a n an- 
alytic solution exist s for this model (originally sta ted bv iLiislI 
Jl952l) and later bv iLvnden-Bell & Pringlel dl974t) L it is fre- 
quently used as a test problem f or numerical codes dev eloped 
to sim ulate accretion discs (e.g.. lSpeith & Rifferll d 1 999h : iKlevI 
Jl999h V 

Usually, modelling the physical shear viscosity by imple- 
menting the viscous stress tensor according to Eq. (1161 is suf- 
ficient to simulate the evolution of the viscous ring. This can 
be seen in Fig. IA.1I where the results of an SPH simulation 
with only physical viscosity, i.e. without any artificial viscosity, 
is shown. Plotted is the azimuthally averaged surface density 
over radius (the radial coordinate is normalized to the radius 
of the initial peak, Rq, and the surface density is normalized 
to the total mass of the accretion ring, M/R^). For the kine- 
matic viscosity we chose v - vq - 10 15 cm 2 /s, the value 
also used in the simulations of the protoplanetary disc. The 
solid line represents the averaged SPH results after a simu- 
lation of 3 x 10 5 sec, the dashed-dotted line the results after 
6 x 10 s sec, while the dashed line indicate the initial distri- 
bution, and the dotted lines denote the analytic solutions, re- 
spectively. The overall time evolution of the ring is reproduced 
very well. The appearing oscillations and deviations from the 
analytic solutions are mainly caused by a secular viscous insta- 
bility inherent to the ring model ( Speith & Kiev 2003). 

However, applying only physical viscosity according to 
Eq. J 16b is insufficient to simulate the dynamical evolution of 
a protoplanet in a disc, because interpenetration of SPH par- 
ticles occurs such that the spiral density wave pattern in the 
disc cannot be resolved. This problem is reliably prevented by 
the new artificial bulk viscosity presented in Eq. d!9i in com- 
bination with the XSPH-scheme of Eq. 11221 . To test for any 
spurious numerical shear viscosity that may be introduced by 
the usage of these algorithms, a simulation similar to that pre- 
sented in Fig. lA.ll was performed with additional artificial bulk 
viscosity, where we set / = 0.5 and XSPH is switched on with 
x - 0.5. The results are shown in Fig. IA.2l where the quantities 
equivalent to those in Fig. lA.ll are plotted. As can be seen, the 
global evolution of the viscous ring, which is very sensitive to 
the effective shear viscosity in the simulation, seems not to be 
affected at all by the new artificial viscosity approach. 

This can be tested further by a simulation where only arti- 
ficial bulk viscosity is used. The results are shown in Fig. IA.3l 
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Fig.A.2. Evolution of the azimuthally averaged surface den- 
sity of the viscous ring at three different times, additionally us- 
ing artificial bulk viscosity and XSPH together with the physi- 
cal viscosity. Dotted lines denote the analytic solutions. 



Any evolution of the viscous ring away from the initial distribu- 
tion, which is indicated by the dashed line, has to be due to in- 
herent numerical shear in the simulation. It can clearly be seen 
that the intrinsic effective shear of the artificial bulk viscosity 
approach is very small. A more quantitative analysis gives a 
value of less than 9 % of vo, which can be neglected in all our 
simulations of the protoplanetary disc. Note that the numerical 
effective shear viscosity is not constant throughout the ring but 
seems to be higher at the edges, especially at the inner edge, 
than in the centre. 

It has turned out that it is necessary to use the XSPH- 
algorithm in addition to the artificial bulk viscosity to prevent 
mutual particle penetration more effectively. This can be seen 
in Fig. IA.4I where two different simulations of the reference 
model of the protoplanetary disc are shown, in the upper panel 
with XSPH switched on and in the lower panel with XSPH 
switched off. The surface density is plotted in greyscale, simi- 



lar to the result presented in Fig. [2 but only for the inner part 
of the disc and in Cartesian coordinates. Both plots show ba- 
sically the same features, but in the case with XSPH the struc- 
tures are more distinct and better resolved, while in the case 
without XSPH the distribution is more diffuse. 

For the sake of completeness it should be mentioned, that 
a frequently used very efficient different approach to prevent 
particle interpenetratio n is the u s ual a rtificial viscosity intro- 
duced by Monashan & Gi n gold! dl983h . and there especially 
the quadratic /?-term. However, this artificial viscosity suffers 
from a high level of spurious intrinsic shear. This becomes ob- 
vious in the simulation of the viscous ring shown in Fig. IA.5I 
where only this artificial viscosity is used with a rather small 
value of p — 1 . The effective physical shear is already of the 
order of vq. 



References 

Artymowicz, P. & Lubow, S. H. 1996, ApJ, 467, L77 
Beckwith, S. V. W., Henning, T, & Nakagawa, Y. 2000, 

Protostars and Planets IV, 533 
Beckwith, S. V. W. & Sargent, A. I. 1996, Nature, 383, 139 
Benz, W. 1990, in Numerical Modelling of Nonlinear Stellar 

Pulsations Problems and Prospects, 269 
Boss, A. P. 2001, ApJ, 563, 367 

Boss, A. P. 2003, in Lunar and Planetary Institute Conference 
Abstracts, 1075 

Bryden, G., Chen, X., Lin, D., Nelson, R. P., & Papaloizou, 

J. C. 1999, ApJ, 514, 344 
Bryden, G., Lin, D. N. C, & Ida, S. 2000, ApJ, 544, 481 
Bubeck, T, Hipp, M., Hiittemann, S., et al. 1999, Parallel SPH 
on Cray T3E and NEC SX-4 using DTS, ed. E. Krause & 
W. Jager (Berlin, Heidelberg, New York: Springer) 
D'Angelo, G., Henning, T, & Kley, W. 2002, A&A, 385, 647 
D'Angelo, G., Kley, W., & Henning, T. 2003, ApJ, 586, 540 
Eggleton, P. P. 1983, ApJ, 268, 368 



Christoph Schafer et al,: Simulations of planet-disc interactions using SPH 



11 



with XSPH 




-4 -2 2 4 
^-coordinate [AU] 



without XSPH 




Fig. A.4. Greyscale plots of the surface density of the inner part 
of the protoplanetary disc in the case of the reference model 
similar to Fig. but in Cartesian coordinates, with (upper 
panel) and without (lower panel) using the XSPH-algorithm. 
Note that the planet cannot be seen on these plots, since it is 
located at a radial distance of 5.2 AU. 



Flebbe, O., Muenzel, S., Herold, H., Riffert, H., & Ruder, H. 

1994, ApJ, 431,754 
Gingold, R. A. & Monaghan, J. J. 1977, MNRAS, 181, 375 
Goldreich, P. & Tremaine, S. 1979, ApJ, 233, 857 
— . 1980, ApJ, 241, 425 

Ivanov, P. B., Papaloizou, J. C. B., & Polnarev, A. G. 1999, 

MNRAS, 307, 79 
Kallenbach, R., Benz, W., & Lugmair, G. W. 2000, in From 

Dust to Terrestrial Planets, 1 
Kley, W. 1999, MNRAS, 303, 696 
— . 2000, MNRAS, 313, L47 

Kley, W., Peitz, J., & Bryden, G. 2003, astro-ph/03 1(5321, 
A&A, in press 

Korycansky, D. G. & Papaloizou, J. C. B. 1996, ApJS, 105, 181 
Lee, M. H. & Peale, S. J. 2002, ApJ, 567, 596 
Lin, D. N. C. & Papaloizou, J. 1985, in Protostars and Planets 
II, 981-1072 

Lin, D. N. C. & Papaloizou, J. C. B. 1993, in Protostars and 

Planets III, 749-835 
Lin, D. N. C, Papaloizou, J. C. B., Terquem, C, Bryden, G, & 

Ida, S. 2000, Protostars and Planets IV, 1 1 1 1 
Lubow, S., Seibert, M., & Artymowicz, P. 1999, ApJ, 526, 1001 
Lucy, L. B. 1977, AJ, 82, 1013 



CO 




radius [R ] 

Fig. A.5. Evolution of the azimuthally averaged surface den- 
sity of the viscous ring at three different times, applying only 
the usual artificial viscosity. The dotted line represents the ini- 
tial analytic solution. 



Lufkin, G, Quinn, T, Wadsley, J., Stadel, J., & Governato, F. 

2003, ArXiv Astrophysics e-prints, 5546 
Lust, R. 1952, Z. Naturforschg., 7a, 87 
Lynden-Bell, D. & Pringle, J. E. 1974, MNRAS, 168, 603 
Marcy, G. W., Cochran, W. D., & Mayor, M. 2000, Protostars 

and Planets IV, 1285 
Mayer, L., Quinn, T, Wadsley, J., & Stadel, J. 2002, Science, 

298, 1756 

Mayor, M. & Queloz, D. 1995, Nature, 378, 355 
Monaghan, J. J. 1989, Journal of Computational Physics, 82, 1 
Monaghan, J. J. 1992, Annual Review Astronomy Astro- 
physics, 30, 543 
Monaghan, J. J. & Gingold, R. A. 1983, Journal of Computa- 
tional Physics, 52, 374 
Nelson, A. F, Benz, W., Adams, F. C, & Arnett, D. 1998, ApJ, 
502, 342 

Nelson, R. P. & Papaloizou, J. C. B. 2002, MNRAS, 333, L26 
Nelson, R. P., Papaloizou, J. C. B., Masset, F, & Kley, W. 2000, 

MNRAS, 318, 18 
Ohtsuki, K, Stewart, G. R., & Ida, S. 2002, Icarus, Volume 

155, Issue 2, pp. 436-453 (2002)., 155, 436 
Papaloizou, J. & Lin, D. N. C. 1984, ApJ, 285, 818 
Papaloizou, J. C. B., Nelson, R. P., & Masset, F. 2001, A&A, 

366, 263 

Papaloizou, J. C. B., Terquem, C, & Nelson, R. P. 1999, in 
Astrophysical Discs - an EC summer school, Conference 
Series, ed. J. A. Sellwood & J. Goodman, Vol. 160 (Astro- 
nomical Society of the Pacific), 186 

Perryman, M. A. C. 2000, Reports of Progress in Physics, 63, 
1209 

Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, 
Icarus, 124, 62 

Sekiya, M., Miyama, S. M., & Hayashi, C. 1987, Earth Moon 

and Planets, 39, 1 
Shakura, N. & Sunyaev, R. 1973, A&A, 24, 337 
Speith, R. & Kley, W. 2003, A&A, 339, 395 
Speith, R. & Riffert, H. 1999, JCAM, 109, 231 



12 



Christoph Schafer et al.: Simulations of planet-disc interactions using SPH 



Syer, D. & Clarke, C. J. 1995, MNRAS, 277, 758 

Terquem, C, Papaloizou, J. C. B., & Nelson, R. P. 2000, in 

From Dust to terrestrial planets, ed. W. Benz, R. Kallenbach, 

G. Lugmair, & F. Podosek 
Ward, W. R. 1997, Icarus, 126, 261 

Lufkin, G., Quinn, T., Wadsley, J., Stadel, J., & Governato, F. 

2004, MNRAS, 347, 421 
Kley, W., Peitz, J., & Bryden, G. 2004, A&A, 414, 735 



